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We reveal a new type of supercritical behavior in gapped graphene with two oppositely charged 
impurities by studying the two-dimensional Dirac equation for quasiparticles with the Coulomb 
potential regularized at small distances accounting the lattice effects. By utilizing the variational 
Galerkin-Kantorovich method, we show that for supercritical electric dipole the wave function 
of the electron bound state changes its localization from the negatively charged impurity to the 
positively charged one as the distance between the impurities changes. Such a migration of the 
wave function corresponds to the electron and hole spontaneously created from the vacuum in bound 
states screening the positively and negatively charged impurities of the supercritical electric dipole, 
respectively. We generalize our results to a particle-hole asymmetric case, where the charges of 
impurities differ in signs and absolute values and demonstrate that the necessary energetic condition 
for the supercriticality of novel type to occur is that the energy levels of single positively and 
negatively charged impurities traverse together the energy distance separating the upper and lower 
continua. The robustness of the supercriticality of novel type is confirmed by the study of an exactly 
solvable ID problem of the Dirac equation with the square well and barrier potential modeling an 
electric dipole potential. 

PACS numbers: 81.05.ue, 73.22.Pr 


I. INTRODUCTION 

It is well known that the Dirac Hamiltonian for the electron in the Coulomb field of a point charge Ze is not self- 
adjoint for Z > 137 and the energy of the lowest 15'i/2 bound state E = my/l — Z'^a^ becomes imaginary testifying 
the fall into the center (atomic collapse) phenomenon^P— . Pomeranchuk and Smorodinsky showed^ that this problem 
disappears if a finite size of nuclei is taken into account. Then a physically acceptable solution exists for larger values 
of the charge and its energy is real. Still the lowest energy electron bound state dives into the lower continuum for 
Z > 170 leading to the spontaneous creation of electron-positron pairs with the electrons screening the positively 
charged nucleus and the positrons emitted to infinity^i^. Since supercritically charged nuclei are not encountered in 
nature, this phenomenon was never observed in quantum electrodynamics. 

It is well known that quasiparticles in graphene are described by the two-dimensional Dirac equation and their 
interaction with the electromagnetic field is characterized by the large effective coupling constant ag = /{hvp) ~ 2.2, 
vp ~ c/300 being the Fermi velocity (c is the velocity of light). Therefore, the value of the critical charge in graphene 
dramatically decreases and equals Zc ~ 1/2^“—. We would like to mention also that the supercritical Goulomb center 
instabiliW is closely related to the excitonic instability in graphene in the strong coupling regime Og > Oc ^ 1 (see 
Refs. [3-0) and possible gap opening, which may transform graphene into an insulator—“i^. 

One would think that the supercritical Coulomb center instability due to the large value of the coupling constant 
should be easily observed in graphene. However, it is difficult in practice to produce highly charged impurities. 
In addition, the external charge in a realistic experimental set-up should be smeared over a finite region of the 
graphene plane because, otherwise, the Dirac equation is no longer applicable and other nearest cr-bands should be 
included in the analysis^. Therefore, the experimental observation of the supercritical instability in graphene was not 
demonstrated until recently. A clever means to solve this problem was recently proposed and realized experimentally 
in Ref. . By creating artificial nuclei in a certain region of graphene fabricated through the deposition of charged 
calcium dimers on graphene with the tip of a scanning tunneling microscope, the supercritical regime was reached 
and the resonances corresponding to the atomic collapse states were observed. 

By making use of the density functional theory and an improved Huckel model, the supercritical instability for 
Ga dimers on graphene was theoretically studied in Ref.[l^. An “atomic-collapse” state in graphene was found for 
fewer absorbed Ga dimers than in the experiment, possibly due to the different spacing between dimers and the 
dielectric screening by a boron nitride substrate. In the continuum model, the study of the supercritical instability 
of one Coulomb center in gapped graphene was extended by usi^ to the case of the simplest cluster of two equally 
charged impurities when the charges of impurities are subcritical, whereas their total charge exceeds a critical one. 
We determined the critical distance between the impurities separating the supercritical and subcritical regimes as a 
function of charges of impurities and a gap. 
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An interesting electric dipole problem in gapped graphene with two oppositely charged impurities was recently 
considered in Refs. [Ml (the 3D Dirac equation with the electric dipole potential was also studied some time ago 
in Ref. El)- It was shown that the point electric dipole potential accomodates towers of infinitely many bound 
states exhibiting a universal Efimov-like scaling hierarchy and at least one infinite tower of bound states exists for an 
arbitrary dipole strength. Notice that the Schrddinger equation in two dimensions for the electron in the field of an 
electric dipole also admits a bound state for any dipole strength^S unlike the three-dimensional case where a bound 
state exists only when the dipole moment exceeds a certain critical value (see, e.g., a discussion including historical 
one in Ref. [2l|). By combining analytical and numerical methods, Ref. [l3| found that the bound states do not dive 
into the lower continuum because the positive and negative energy levels first approach each other and then go away. 
Actually this behavior is typical for an avoided crossing^^, which forbids level crossing for two states with the same 
quantum numbers. Since the bound states do not dive into the lower continuum, the authors of Ref. (a concluded 
that supercriticality is unlikely to occur in the electric dipole problem in graphene. 

We reconsidered the problem of supercriticality in our recent paper [23| for the case of two oppositely charged 
impurities situated at finite distance (finite electric dipole). By using the linear combination of atomic orbitals (LCAO) 
and variational Galerkin-Kantorovich (GK) methods, we showed that for sufficiently large charges of impurities the 
wave function of the highest energy occupied bound state changes its localization from the negatively charged impurity 
to the positive one as the distance between the impurities changes (both methods gave similar results). The necessary 
condition for the instability to occur is the crossing of the electron energy levels in the field of single positively and 
negatively charged impurities. This migration of the electron wave function of the supercritical electric dipole is a 
generalization of the familiar phenomenon of the atomic collapse of a single charged impurity with holes emitted to 
infinity to the case where both electrons and holes are spontaneously created from the vacuum in bound states with 
two oppositely charged impurities thus partially screening them. 

In this paper, we extend by using the GK method the study of the electric dipole problem in gapped graphene 
performed in Ref. . We apply this method also to the more general case of two oppositely charged impurities whose 
charges are not equal by modulus. In this case, we find that the wave function of the highest energy occupied bound 
state changes its localization only if the bound state levels in the corresponding one Goulomb center problems traverse 
together the energy distance 2 A (A is a gap). Since the LG AO and variational GK methods are approximate ones 
because in practice one can take into consideration only a finite number of trial functions, it is important to check the 
robustness and validity of the supercriticality of novel type connected with the migration of the wave function of the 
electron bound state. For this, we study in this paper an exactly solvable ID model of the Dirac equation with the 
square well and barrier potential modeling an electric dipole potential. In addition, this model allows one to consider 
large constituent charges in a dipole when oscillations of energy levels and connected with them the changes of the 
localization of the electron wave function take place. 

The paper is organized as follows. The supercritical instability in an exactly solvable ID problem is investigated 
in Secini In Sec. IIIIl we consider the symmetry properties of the Dirac equation for quasiparticles in the field of an 
electric dipole in graphene. In Sec lIVl we solve the Dirac equation and study the supercritical instability in graphene 
in the electric dipole potential by means of the numerical variational Galerkin-Kantorovich method. An asymmetric 
problem of two impurities with different by modulus charges of opposite sign is considered in SeclV) The results are 
discussed in Gonclusion. Oscillations of energy levels in the exactly solvable ID problem are considered for sufficiently 
large strength of an electric-dipole-like potential in Appendix The system of differential equations in the variational 
Galerkin-Kantorovich method is written down in Appendix [BJ 


II. ONE-DIMENSIONAL MODEL WITH AN ELECTRIC-DIPOLE-LIKE POTENTIAL 

In this section, we consider the supercriticality of novel type and the effect of relocalization of the wave function of 
the highest energy occupied bound electron state in an exactly solvable onc-dimensional model of the Dirac equation 
with the square potential well and barrier modeling an electric dipole potential^. This model can be applied also for 
the description of edge states in graphene in the presence of a dipole on a boundary. The Hamiltonian of this problem 
reads 


H{A,Vo) =-ihvpCTxdx + Aa^+ V{x), (1) 

where and cr^ are the Pauli matrices, A is a gap, and an “electric-dipole-like” potential V(x) is defined by the 
equation 

V(x) = —VqO {d—\x + a|) -I- VqO {d — \x — a|) (2) 

and is plotted in FiglT] Here 2d is the width of the well and barrier, and 2a is the distance between their centers. In 
what follows, we will use dimensionless quantities e = E/A, vq = Vq/A, x —> xRa = xhvp/A, a —> ahvp/A, and 




d —>■ dhyp/^ (note that hyp/^ is the Compton wavelength). Hamiltonian ([T} with potential (l2|) has a particle-hole 



FIG. 1: (Color online) An “electric-dipole-like” potential of the ID model. 

symmetry expressed by HH(A, Vb)H+ = —H{A, Vq), where the unitary operator H = axTZx, with the operator TZx of 
reflection x —>■ —x, satisfies = 1. It follows then that all solutions of the Dirac equation come in pairs with ±e. 

We will see below that the behavior of the energy levels and wave functions in the model with the “electric-dipole- 
like” potential depends on the sign of energy of the electron bound states in the single potential well and barrier 
problems. For the potential v{x) = —VgO {d — |a;|), we plot in Figl2]the dependence of the energy levels of electron 
bound states on |no| for the potential well and barrier at the fixed value d = 0.25. The energy levels in the potential 
well with uq > 0 are plotted by solid blue lines and the corresponding levels for the potential barrier of the same 
strength but opposite sign potential uq < 0 are shown by red dashed lines. As is seen, bound states appear for an 
arbitrary small |no|. As |no| increases, the energy levels cross e = 0 and then enter the continua for certain critical 
strengths (for example, |no| Ri 7.5 for the first levels). The energies of the first bound state levels change their sign at 
I'col ~ 4. 
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FIG. 2: (Color online) The dependence of energy of bound states in the potential well/barrier on the strength of the potential 
|tio| at the fixed value d — 0.25: blue solid lines correspond to the potential well (no > 0) and red dashed lines to the barrier 
(no < 0). 

We will solve the Dirac equation with the electric-dipole-like potential in each of five regions defined by Eq. ([2l) and 
then match solutions at the points where the potential jumps by using the condition of continuity of each component 
of the spinor. Then, for the component (j) of the two-component spinor I'F) = {(j>, x)^, we find the equation 

(j)" + {{v - e)2 - I) 0 = 0, 

and X is related to (f) through y = Wxip/(v — e — 1). We have the following solutions: 

I) a; < —a — d 


( 3 ) 
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2) —a — d < X < —a + d 


3) —a + d<x<a — d 


4) a — d<x<a + d 


5) X > a + d 


( 1)11 (x) 
XII(x) 


C 2 sin(fcix) + C 3 cos(fcix), 
ki 

—i —q-— {C 2 cos(fcix) — C 3 sin(fcia;)); 


cj)ui{x) = , 

Xiii(x) = (C'4e"“" - C 5 e-««") ; 

(t)iv{x) = Ce sin(fc2a;) + Cy cos(fc2a:), 

^2 

Xiv(x) = i -- (Ce cos(fc 2 a;) — C 7 sm{k2x)) ; 

vq - e - 1 

cj^vix) = Cse-^°\ 

XV {x) = 


( 5 ) 


( 6 ) 


( 7 ) 


( 8 ) 


Here we defined Kq = \/l — e^, ki = ■>/(vg + e)^ — 1, and /c 2 = (vg — e)^ — 1. By matching the solutions 

4>(xi + 0) = (^(xi - 0), x(xi + 0) = x(xi - 0), i = l,4 (9) 

at the four points 

xi = —(a + d), X 2 = —ia — d), X 3 = a — d, X 4 = a + d, (10) 

we find the following system of equations for coefficients Ck- 

A,kCk = 0 , ( 11 ) 

where the coefficients Aik are given by Ea. dAll) in Appendix A. The secular equation det A = 0 determines the bound 
state energy levels of the problem. For a < d, the energy levels are obtained by the interchange a d. The coefficients 
Ck determined from Eq. (HU together with a normalization condition specify completely the wave functions. 

For d = 0.25, we plot the bound state energy levels in Figl^for two values z;o = 3 and vq = 6. 
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FIG. 3: (Color online) The bound state energy levels of the ID problem with the electric-dipole-like potential for vg = 3 (red 
dashed lines) and uo = 6 (blue thick lines) obtained for d = 0.25. 

Figl3] implies that the energy levels for the smaller value uq = 3 monotonously depend on the distance between the 
well and barrier. We plot in FigHlthe square modulus of the wave function of the negative energy bound state for 
two values of the distance between the centers of the well and barrier 2a = 0.3 and 2a = 1.5 whose potentials are 
schematically shown by filled green regions. Since the well and barrier overlap for 2a = 0.3, their widths are reduced 
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in this case. For 2a = 1.5, the well and barrier are well separated and, obviously, the wave function does not change 
its localization on the barrier. Note that for Uq = 3 the energy levels in the single potential well problem do not cross 
the zero energy value (see Fig. E]). 

For the larger value uq = 6, according to FigEl the level repulsion is observed. We plot in FigOthe square modulus 
of the wave function of the negative energy bound state for four values of the distance 2a between the centers of the 
well and barrier. For 2a = 0.14, the wave function is localized on the barrier whose width is much reduced due to the 
strong overlap with the well. For larger value of the distance between the centers of the well and barrier 2a = 0.35, 
the wave function is localized both on the barrier and well. As the distance between the centers increases further, the 
wave function migrates to the well. It is crucial that the case uq = 6 corresponds to the negative sign of the energy of 
the first bound state in the potential well problem (see Fig. [2]). Note that a similar phenomenon of the change of the 
localization of wave function takes place in the fission of quarkonium consisting of heavy quark and antiquar k^i^^ . 




FIG. 4: (Color online) The square modulus of the wave function of the negative energy bound state for vq = 3, d = 0.25, and 
two values of the distance between the centers of the well and barrier: (a) 2a = 0.3, and (b) 2a = 1.5. The potentials of the 
well and barrier are schematically plotted as filled green regions. 
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FIG. 5: (Color online) The square modulus of the wave function of the negative energy bound state for vq = Q, d = 0.25, and 
four values of the distance between the centers of the well and barrier: (a) 2a = 0.14, (b) 2a = 0.35, (c) 2a = 0.7, and (d) 
2a = 1.6. The potentials of the well and barrier are schematically plotted as filled green regions. 

For larger values of uq when several energy levels cross the zero energy e = 0 in the potential well problem in FigEJ 
we observe oscillations in the behavior of energy levels as well as oscillations in the localization of the wave function 
on the well and barrier (see. Appendix . 

Thus, the exact solutions of the ID Dirac equation with the electric-dipole-like potential unambiguously confirm 
the conclusion made in Ref. that the supercritical instability in the presence of both attractive and repulsive 
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potentials is connected with the change of the localization of the wave function of the highest occupied electron bound 
state. 

The local density of states LDOS(e,a,a;) = J2k where k is the set of all quantum numbers, is 

a physical quantity that can probe the migration of the wave function and be directly measured in an experiment. 
Therefore, we present the video in Supplemental Material where the LDOS in the ID model with an electric- 
dipole-like potential considered in this section is plotted as a changes for the energy close to the boundary of the 
lower continuum e = —1.1 and the three values of vq such that the energy of the lowest electron bound state in 
the potential well problem eo equals 0.35, —0.085, and —0.88. Obviously, the first case is subcritical because eo is 
positive, therefore, the peak of the LDOS remains localized on the barrier as the distance 2a between the centers of 
the well and barrier changes from small to large values (note that since the peak is localized on the barrier, it moves 
to the right as 2a increases). The two other cases are supercritical and the migration of the peak of the LDOS from 
the barrier to well is clearly seen as a changes. Since the last case eo = —0.88 is strongly supercritical, the LDOS 
peak takes much larger values in this case compared to those in the weakly supercritical case eo = —0.085 when the 
lowest energy bound state level in the square well potential is only slightly below the zero energy. Thus, we conclude 
that measuring the LDOS in the continua makes it possible to demonstrate the supercriticality of novel type and the 
migration of the electron wave function in the electric dipole problem. 


III. DIRAC EQUATION 

Let us consider now the supercritical instability in the electric dipole problem in gapped graphene. The Dirac 
Hamiltonian in 2 -|- 1 dimensions which describes the quasiparticle states in the vicinity of the K± points of graphene 
in the field of two oppositely charged impurities reads (we set fi = 1) 

H{A,e) = vpcrp + - eV{r), (12) 

where — e < 0 is the electron charge, p = —iVr is the two-dimensional canonical momentum, ai are the Pauli 
matrices, and A is a quasiparticle gap. The latter can be opened in graphene in various ways, e.g., due to finite- 
size effects in graphene nanoribbons^^ or by depositing graphene on a substrat o^^i^° . Hamiltonian (HI acts on two 
component spinor which carries the valley (^ = ±) and spin (s = ±) indices and we use the standard convention: 
'i'+s = {'4’A,iJB)K+s, whereas = {'ipB,'4’A)K-s, and A,B refer to two sublattices of hexagonal graphene lattice. 
We regularize the Coulomb potential of each impurity by rp, which is of the order of the graphene lattice spacing. 
Then the regularized interaction dipole potential for charged impurities AQ, Q = 2'e, situated in the (x, y) plane at 
(±i?/2,0) is given by 


vi^) = -( ^ ^ 

K \ 1 /( 3 ; -K i?/2)2 +y2 +r‘^ \/{.X- i?/2)2 + ^2 

where k is the dielectric constant. For the sake of definiteness, we will consider the electrons only in the Kj^ valley 
(the Dirac equation for the electrons in the K- valley is obtained by replacing A with —A). Since the interaction 
potential (1131) does not depend on spin, we will omit the spin index s in wave functions in what follows. The main 
difficulty in solving the Dirac equation for the electron in the electric dipole potential is that variables in this problem 
are not separable in any known orthogonal coordinate system. Therefore, we will utilize in this paper the numerical 
variational Galerkin-Kantorovich method. 

Let us discuss the discrete symmetries of Hamiltonian (USD at fixed valley and spin with potential (1131) . The 
parity transformation changes A to —A and also reverses the sign of the potential, or equivalently in the case under 
consideration, it reverses the charge of the particle, e —>■ —e: 

t7piJ(A,e)t/-i=iL(-A,-e), (14) 

where Up = UylZx is a unitary operator with the operator TZx of reflection {x,y) -A {—x,y). Hence the wave function 
'l'p(r) = [/p4/(r) = ay^{TZxr) describes states with the same energy but the opposite sign of the charge and gap. 

The charge conjugation operator Uc = <JxK, where K is the complex conjugation, interchanges the Hamiltonians 
H{A,e) and H{A,—e) 



UcH{A, e)U-^ = -HiA, -e). (15) 

Therefore, if the wave function 4'(r) is a solution of the stationary Dirac equation iJ(A, e)4'(r) = i?4'(r), then the 
charge conjugated wave function 4'c(r) = 17c4'(r) = ax'i’*{x) is an eigenfunction of H{A, —e) but with the eigenvalue 
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—E. The time-reversed wave function = C/T'It(r) = crydt*(r) satisfies the equation iJ(—A, e)4'T(i’) = i?4'T(r) 

because 


UtH{A, e)UT^ = H{-A, e). (16) 

This reflects the well known fact that a gap at fixed valley and spin in graphene breaks the time reversal symmetry. This 
symmetry is effectively preserved in the charge density wave and quantum spin Hall states^, which are characterized 
by a gap of opposite sign for the other valley and spin, respectively. 

Hamiltonian (TT^ with potential (1131) has an intrinsic particle-hole symmetry expressed by Hi7(A, e)H+ = —H{A, e), 
where the unitary operator O = axE-x satisfies = 1 (note that O is the same operator as in the ID model in Secllll). 
It follows then that an eigenstate E{x,y) with energy E has a partner '^_E{x,y) = Hd>£;(a;,?/) = ax'^Ei—x^y) with 
energy —E, hence, all solutions of the Dirac equation come in pairs with ±E. In fact, the operator D is nothing else 
as the combination of three symmetry operations C, P, and T: D = UcUpUT- The intrinsic particle-hole symmetry of 
the electric dipole problem will play a prominent role in our analysis below. 

The Dirac Hamiltonian with the electric dipole potential commutes also with the operator U = azKTZy, where 
TZy is the operator of reflection y —>■ —y. The equality = 1 implies that wave functions are split into two classes 
HldtA) = AldtA), where A = ±1. Since the operator U is antilinear, the function |d>_) related to the function |'I'+) 
by means of the phase transformation is an eigenfunction with A = — 1. Therefore, there no need 

to consider functions |'I'_). Hence the components of wave functions (r|d)'+) = '!/i(r) = (<)>, x)^ satisfy the following 
constraint conditions consistent with the Dirac equation: 

<('*(-2/) = <('(2/), /lyN 

x*{-y) =-xiy)- 


It is convenient to work with dimensionless quantities h = ^ and e = In addition, we assume in what follows 
that all coordinates and distances are dimensionless and are defined in units of i?A = We introduce also 

dimensionless coupling constant C = Then the Dirac equation reduces to the following system of two coupled 

ordinary differential equations of the first order: 

( -i{dx + idy)(l) + {v - e - l)x = 0, 

\ -ildx - idy)x+ lv - e + 1)^ = 0. ^ ' 

We will find numerical solutions of Eas. (fTSl) in the next section by using the variational GK method in the class of 
wave functions with A = 1. Still it is instructive to begin our analysis with the Dirac equation for the electron in 
graphene with one positively charged impurity 




(19) 


where 


hp — T ^ydy') ~\~ (Tz 


c 


( 20 ) 


The Hamiltonian for the electron in the field of negatively charged impurity is obtained from the Hamiltonian hp 
by the change of the sign of the last term in hp. 

The wave function of the state with the total angular momentum j in the polar coordinates (r, 9) has the form 


T = 


e*(^-i/2)V(r) \ 

g{r) J 


( 21 ) 


For the lowest energy electron bound state with j = 1/2, the Dirac equation takes the form 


/ /'= (l + e-i^(r))g, 

1 5 ' + 7 = (1 -i + v{r))f. 


( 22 ) 


To find the boundary conditions for the functions / and g at the origin, we investigate the asymptotic behavior at 
r = 0. The solution of the approximate equation 


/" + ^ + \/ = 0 


( 23 ) 





regular at the origin is given by / = Jq (C^/'^’o); where Jo{z) is the Bessel function. Therefore, the numerical solution 
of the Dirac equation should satisfy the boundary conditions /(O) = 1 and (?(0) = 0. 

We determine numerically the energy levels by using the shooting method and requiring that the wave functions 
decrease at infinity. We normalize the wave functions according to the condition f = 1. The energy of the 

lowest (highest) electron bound state with j = 1/2 in the regularized Coulomb potential with the charge +Q (—Q) 
is plotted in Figl 6 ] for several values of ro as a function of The levels which descend from the upper continuum 
correspond to the positively charged impurity +Q (these levels agree with the corresponding results in Ref.Q, see 
Fig.4 there), while those which rise from the lower continuum and grow with ( correspond to the negatively charged 
impurity —Q. These results also reproduce qualitatively the behavior seen directly at the tight-binding level on a 
honeycomb lattice f^]. 

For nonregularized Coulomb potential of positive charge, the energy of the lowest bound state is always posi¬ 
tive. It reaches the value e = 0 for ( = 1/2 and becomes purely imaginary for (> 1/2 (the “fall into the center 
phenomenon”— “—). For regularized Coulomb potential, the energy of the lowest bound state crosses e = 0 and ap¬ 
proaches the negative-energy continuum for a certain critical charge. For example, for ro = 0.05i?A, the critical charge 
() Ri 1 corresponds to the lowest electron bound state diving into the lower continuum. 



FIG. 6: (Color online) The energy of the electron bound state with j = 1/2 in the regularized Coulomb potential as a function 
of ( for different valnes of the regularization parameter: ro = 0 (green dash-dotted lines), ro = 0.01i?A (blue dashed lines), 
ro = 0.05i?A (red solid lines). The levels which descend from the upper continuum correspond to the positive charge +Q 
problem while those which rise from the lower continuum and grow with ( correspond to the negative charge —Q problem. 

The electron levels in the field of negatively charged center described by the Hamiltonian hn are obtained by the 
reflection e —>■ —e because the operator of charge conjugation 14 = interchanges the hp and hn Hamiltonians, 
UchpU^ = —hn- Therefore, the energy levels of the Hamiltonians hp and intersect at e = 0. We will see that 
the corresponding critical value of a coupling constant 4 pl^-ys a crucial role in our analysis of solutions in the 
electric dipole potential because the behavior of the corresponding energy levels dramatically changes depending on 
whether ^ or (^ > 4 similarly to the analysis of the ID model in the previous section. For the chosen values of the 
regularization parameter rg in FiglHl the critical coupling constant Cc — 0.6 {tq = O.OIRa) and Cc — 0.7 (r-g = O.OSRa)- 
Finally, the Hamiltonian for quasiparticles in graphene with two oppositely charged impurities in dimensionless 
units has the form 

h =-i{aA + (Jydy) + + C \ - , ^ + / ■} 9 I ^ (24) 

V Vrl+r^o) 

where Vp^n = \/{x ± 7?/2)^ -I- In order to find eigenstates of this Hamiltonian, we will apply the variational 
Galerkin-Kantorovich method in the next section. 


IV. VARIATIONAL METHOD 


The choice of trial functions is a crucial element in any variational method. We choose the following trial functions 
in the GK method which belong to the class A = 1 (see the discussion above Eg. (1171) 1 and satisfy the asymptotic at 
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large distance: 


(j){x,y) = e fl/ 2 )^+y 2 +rg | ^ f 2 k{x)y^'‘ +i J2 f 2 k-i{.x)y'^'‘ A , 

[fe=0 fc=l J 

X{.x,y) = I ^ g2k{x)y'^^ +i Y, g2k-i{x)y'^^~^'\ ■ 

[fc=0 fc=l J 


( 25 ) 


According to the Galerkin-Kantorovich method, the wave functions are substituted into the initial equation and their 
orthogonality to the residual with respect to the variable y is required. Thus, we obtain the system of ordinary 
differential equations (see Ea. dBlI) in Appendix I bI) for functions fk and gk- 


€ 



FIG. 7: (Color online) The energy levels of bound states obtained in the variational method for ro = O.OSAa and different 
charges of impurities: ^ — 0.6 (green dashed lines) and = 0.85 (blue solid lines) 


To determine how successful the GK method is, we first use this method for solving the Dirac equation for the 


electron in the field of one Coulomb center whose regularized potential is given by V{r) = — 




We should set 


i? = 0 in Eas. dBlI) - dB3l) and take into account that 14(a;) = —CQs{x). The exact spectrum of the one Coulomb center 
problem with ro = 0.05i?A plotted in Figl6](red solid lines) can be compared with the approximate spectra determined 


by the GK method for the various number of terms in the ansatz: TV = 0,7V = 0; TV = l,iV = 0; TV = 0,TV = 1; 
and TV = 1,TV =1. We checked that the approximations “00” and “11” reproduce the exact spectrum better than 
the other approximations. Therefore, we use in what follows these approximations to determine the bound states for 
the electron in the electric dipole potential. 



FIG. 8: (Color online) The square modulus of the wave function of the lower energy level for = 0.6 and various distances 
between the impurities: R = 1.25 (left panel), R = 2.0 (middle panel), and R — 3.0 (right panel). The wave function doesn’t 
change its localization from the negatively charged impurity to the positively charged one as the distance between the impurities 
increases. 

The approximate spectrum of the Dirac equation for the electron in the electric dipole potential with rg = 0.05i?A 
is plotted in Fig. [3for C = 0.6 (green dashed lines) and ( = 0.85 (blue solid lines). For ( = 0.6, the energy of the 
lowest electron bound state in the single positive Goulomb center problem (see FiglH]) is positive. Like in the ID 
model, we find then that the wave function doesn’t change its localization on the impurities (see Fig. El). For C = 0.85, 
the energy of the lowest electron bound state in the single positive Goulomb center problem (see Fig. [6]) is negative 
(for the chosen charge ( = 0.85, the single positive Coulomb center has only one negative energy level). Therefore, the 
wave function changes its localization on the charged impurities (see Fig. E]). Our analysis performed by making use 
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FIG. 9: (Color online) The square modulus of the wave function of the lower energy level for ^ = 0.85 and various distances 
between the impurities: R — 1.0 (left panel), R = 2.55 (middle panel), and R = 3.0 (right panel). The wave function changes its 
localization from the negatively charged impurity to the positively charged one as the distance between the impurities increases. 


of the GK variational method shows that the migration of the wave function of the highest energy occupied electron 
bound state takes place when the charges of impurities are such that the energy levels of the corresponding single 
Coulomb centers with charges ±Q cross. 

For larger more complex behavior of the energy levels compared to that given by the blue solid lines in Fig|7]with 
more cycles of oscillations of energy levels is observed as found in Ref. da and qualitatively similar to that in Fig ll2l 
in the exactly solvable ID model in Appendix]^ It is shown in Appendixthat the oscillations of energy levels in 
the ID model correspond to the oscillations of the localization of the wave function from the potential barrier to well. 


V. ASYMMETRIC CASE 


By using the variational Galerkin-Kantorovich method, we established in the previous section the supercritical 
instability for quasiparticles in graphene in the electric dipole potential. This instability is connected with the change 
of the localization of the wave function of the highest energy occupied electron bound state. We saw that the 
necessary condition for the supercriticality to occur is that the energy of the electron bound state in the field of 
the single positively charged impurity cross zero. Since there is a charge conjugated bound state level in the single 
negatively charged impurity problem, taken together these two levels traverse the energy distance 2A. Recall that the 
supercritical instability in the single Goulomb center problem takes place when the lowest energy electron bound state 
traverses also the distance 2A. This suggests that if we consider the Dirac equation for quasiparticles in graphene 
with two impurities whose charges have opposite sign, however, are not equal by modulus, the supercritical instability 
will take place only if the lowest and highest bound state levels of the corresponding single positively and negatively 
charged impurity problems traverse together the energy distance 2A. We will study this suggestion in this section by 
considering the Dirac equation for quasiparticles in graphene with two impurities whose charges are opposite in sign 
and not equal by modulus. It is obvious that the intrinsic particle-hole symmetry defined by the operator fl in Sec lIIIl 
is no longer present in the case of two oppositely charged impurities with charges not equal by modulus. 

The Hamiltonian of this problem is given by Eg. (1121) with the potential 


V{r) = e 


__ 

y/{x - R/2)'^ + y'^ + 


_ ^^ 

^{x + i?/2)2 + y'^ + r^j 


(26) 


In terms of dimensionless quantities, the Dirac equation for the electrons in the K+ valley has form (fT^ . The Dirac 
Hamiltonian with potential (E51) has the discrete symmetry described by the operator U introduced in Sec lIIIl The 
components of the wave function |'I'+) satisfy conditions (HZ}. In order to solve the Dirac equation with potential 
(|M1) , we use the variational Galerkin-Kantorovich method and utilize the trial wave functions (1251) which satisfy the 
asymptotic at large distance and conditions (HZ}- 

According to the Galerkin-Kantorovich method, we substitute the trial wave functions into the Dirac equation and 
require their orthogonality to the residual with respect to the variable y. We obtain the same system of ordinary 
differential equations as in the symmetric case (IB1|) with the coefficient functions (IB2|) - (IB4|) . The coefficient function 
(IB4I) should be calculated for potential (1^ . The condition that functions fk{x) and gk{x) are finite as x —)► ±oo 
allows us to determine the spectrum. 

In Figlini we plot the bound state energy levels in two cases Ci = 0.5, ^2 = 0.8 (red dashed fines) and Ci = 0.7, 
C 2 = 0.95 (blue solid lines) for tq = 0.05Ra found by using the variational method with N = 0,N =0 (one term in 
the ansatz). Note that the energy levels are not particle-hole symmetric. As i? —>■ 0, the energy levels correspond to 
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FIG. 10: (Color online) The energy levels in the problem of two oppositely charged impurities with ((i = 0.5, (^2 ~ 0.8 (red 
dashed lines) and ((1 = 0.7, ^2 = 0.95 (blue solid lines) for ro = 0.05i?A found in the variational method. 



FIG. 11: (Color online) The square modulus of the wave function of the lower energy level for ((1 = 0.7, (2 = 0.95 and 
various distances between impurities: R = 2.0 (left panel), R = 3.25 (middle panel), and R = 3.75 (right panel). The wave 
function changes its localization from the negatively charged impurity to the positively charged one as the distance between 
the impurities increases. 


the levels of single positive charge impurity with Z 2 — Zi > 0. For R —>■ 00 , the energy levels tend to those related 
with single positive and negative charge impurities. The negative energy levels arise from the lower continuum at 
some critical distance Rc, and this happens only when two impurities become sufficiently separated. 

In the case Ci = 0-5 and (2 = 0.8, the bound state energy levels monotonously depend on the distance between 
impurities and no level repulsion is observed. We checked that the wave function does not change its localization 
similar to the case of two oppositely charged impurities whose charges are equal by modulus. For larger values Ci = 0.7 
and (2 = 0.95, the energies of the bound state energy levels given by the blue solid lines in Fig. (TU] first converge 
and then go away from each other. This suggests that the wave function of the lower energy bound state should 
change its localization. In Fig. 1111 we plot the square modulus of the wave function for the lower energy level at 
three different distances between the impurities. Clearly, the wave function of the highest energy occupied state does 
change its localization from the negatively charged impurity to the positively charged one as the distance between 
the impurities increases. Thus, the study of an asymmetric case in the present section confirms the universality of 
the phenomenon of the change of localization of the wave function and the necessary condition for the supercritical 
instability to occur is that the lowest and highest energy bound states of the corresponding single Coulomb center 
problems traverse together the energy distance 2A. 


VI. CONCLUSION 

Motivated by a recent study of the Dirac equation for quasiparticles in the electric dipole potential in graphene, 
we studied in the continuum model the supercritical instability in the electric dipole problem in gapped graphene. 
Since the variables for the Dirac equation with the electric dipole potential are not separable in any known orthogonal 
coordinate system, we used the numerical variational Galerkin-Kantorovich method with trial functions which have 
the correct asymptotic at large distances. For sufficiently large charges of impurities such that the lowest energy 
electron bound state of the single Coulomb center problem crosses the zero energy E = 0, the positive and negative 
energy bound state levels for the Dirac equation in the electric dipole potential at first converge and then go away 
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(since the levels have the same quantum numbers, they do not cross due to the avoided crossing theorem) tending 
to the energy levels of the corresponding single Coulomb center problems. We found that the wave function of 
the highest energy occupied bound state level changes its localization from the negatively charged impurity to the 
positively charged one at the point of the closest convergence of the bound state levels. This migration of the wave 
function corresponds to a supercriticality of novel type with the spontaneous creation of an electron-hole pair in bound 
electron and hole states screening the negatively and positively charged impurities, respectively. 

We extended our analysis of the electric dipole problem in graphene to an asymmetric case, where the charges 
of impurities are of opposite sign and not equal by modulus. We found in this case that the wave function of the 
highest energy occupied bound state changes its localization only if the bound state levels in the corresponding single 
Coulomb center problems traverse together the energy distance 2 A separating the upper and lower continua. The 
supercritical instability of the electric dipole can be observed experimentally by placing subcritical oppositely charged 
impurities on graphene and then moving with the tip of scanning tunneling microscope one impurity toward the 
other one and afterwards moving the impurities apart again. The supercritical instability takes place if the impurities 
become screened. 

Since the LCAO and variational Galerkin-Kantorovich methods are approximate ones, in order to study the robust¬ 
ness and validity of the supercriticality of novel type we studied an exactly solvable ID model of the Dirac equation 
with a square potential well and barrier modeling an electric dipole potential. Our findings in this exactly solvable 
ID model unambiguously demonstrate the presence of the supercriticality of novel type connected with the migration 
of the wave function of the electron bound state. 
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Appendix A: Oscillations of energy levels in ID model with an electric-dipole-like potential 


In the ID Dirac model 0 with the electric-dipole-like potential ([2]), the energy spectrum is determined from the 
secular equation det A = 0, where the matrix elements are given by 


All 

A 21 

A 32 

A 54 

A 42 

A 45 

Aye 

Age 


-Ays = e A12 = sin [(a-I-(i)fci], A13 = - cos [(a-I-(i)fci] 


Koe 


— (a+d)Ko 


1 + e 

— sin [(a — d)ki ], 
-A 35 

fei cos [(a — d)ki] 


A 22 — — 


ki cos [(a -I- d)ki 


A 23 = — 


ki sin [(a -I- d)ki 


1 + e + VQ 1 + e + VQ 

A33 = cos [(a - d)ki ], A55 = -A34 = 

A56 = - sin [(a - d)k2 ], Agy = - cos [(a - d)k2 ], 


1 -I- e -I- uo 


A 43 = 


fci sin [(a — d)ki 


A64 = 


Koe 


{a—d)Ko 


1 + e 
sin [(a -I- d)k2 ], 
k2 cos [(a -I- d)k2] 
I+ e-Vo 


Aee = — 


1 -I- e -b uo 
k2 cos [(a — d)k2] 


1 + e — Vo 

Ayy = COS [(a -b d)k2] , 


A 44 = Aes = — 
, Aey = 


Koe 


— {a—d)KQ 


l + e 
/c2 sin [(a — d)k2] 


l + e-Vo 


Agy = — 


k2 sin [(a -b d)k2] 
1 + e-Vo 


(Al) 


and the rest of elements equals zero. 

We analysed in Sec|lT]the behavior of energy levels for sufficiently small value Vo- It is instructive to consider what 
happens for larger values of vq. We plot in Fig. [ 12 ] the bound state energy levels in the “electric-dipole-like” potential 
for d = 0.25 and two values of the strength of the potential: uo = 9 (green dashed lines) and vq = 12 (blue solid 
lines). According to Fig. | 2 l for uo = 9 , only one energy level in the potential well problem crosses the zero energy 
e = 0 . This corresponds to the only point of level repulsion (the point of maximal convergence of levels) in Fig. [T 2 | 
(green dashed lines). For vq = 12 , two levels in the potential well cross the zero energy e = 0 (see Fig. | 2 |). This 
corresponds to the two points of level repulsion in Fig. [ 12 ] (blue solid lines) and between these two points there is a 
point of maximal divergence of energy levels. Thus, a qualitatively new phenomenon connected with oscillations of 
energy levels appears for sufficiently large values of Vq. 
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FIG. 12: (Color online) The bound state energy levels of the ID problem for uo = 9 (green dashed lines) and vo = 12 (blue 
solid lines) obtained for d = 0.25. 


In order to demonstrate how the oscillations of energy levels for vq = 12 are reflected in the localization property 
of the wave function of the highest energy bound state, we plot in Fig. [T3]the square modulus of the electron wave 
function of the negative energy bound state at seven values of a shown by red points in Figll2l One can see in 
FigHSl (b) that the electron wave function is localized both on the barrier and well at the first point of the maximal 
convergence of the energy levels. According to Fig ll3l (d), the wave function at the point of the maximal divergence 
of levels between the points of the maximal convergence is localized mainly on the well. At the second point of the 
maximal convergence of the energy levels, the wave function whose square modulus in plotted in Fig ]13l (f) is again 
localized both on the well and barrier. As the distance between the centers 2a increases further, the electron wave 
function becomes localized on the well. 

For still larger values of vq, we observe more oscillations of energy levels which qualitatively resemble those found 
in Ref. El for the electric dipole problem in graphene. 


Appendix B: Differentional equations in the Galerkin—Kantorovich variational method 


For the trial functions fk and gk in Ea. (l25l) . we find the following system of equations {I = Q,N; I =1,N'): 

|Pfe+i - (1 + e)52fe(a;)^ - \/l - - ^sign(a;)^ Qk+if2k{x) + Vk+ig2k{x)'^ - 

n' 

y] {(2A: - l)Pk+i-i - Vl-e^Qk+i] f2k-iix) = 0, 

df 2 k-i{x) e)g 2 k-i{x)^ - - ^sign(a;)^ Qk+i'-if 2 k-i{x) + Rfc+y _i 52 fc-i(a;)| + 


/c=0 

N' 




/c=l 


\ Pk+i'-i 


dx 


+ yi d^Qk+i'] hk{x) =0, 

|pfc+; - (1 - e)/ 2 fc(a;)^ - \/l - - ^sign(x)^ Qk+ig 2 k{x) - Vfc+z/ 2 fe(x)| + 


N' 


+ y^ {(2A:- l)Pk+i-i - \/l - g 2 k-i{x) = 0, 




^ ^Pk+i'-i “ e)/2fe-i(a;)^ - Vl - (^x - ^sign(x)^ Qk+i'-i92k-iix) - _i/2fc-i(a;)| - 


N 

X! £^Qk+l'] 92k{x) = 0, 


(Bl) 


k^O 
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FIG. 13: (Color online) The square modulus of the wave function of the negative energy bound state for vo = 12, d — 0.25, 
and seven values of the distance between the centers of the well and barrier: (a) 2o = 0.08, (b) 2a = 0.144, (c) 2a = 0.19, (d) 
2o = 0.25, (e) 2a — 0.34, (f) 2a — 0.42, (g) 2a = 1.0. The corresponding values of a are marked by red points in Fig. 1121 The 
potentials of the well and barrier are schematically plotted as filled green regions. 


where the coefhcient functions equal 

+ 00 


F4x)= [ e-2vT3?^(|^|-fl/2)^+y2+,2^2,^^ ^ x^+\2s 

J o'* 


+ 00 


Q,(x) = / e-^VT^y/(\-\-m?+y^+rl 


|a;| - RI2) +y'^ + rg 


--dy = ---Ks(ax), 


(B2) 


(B3) 


Vs{x) 


+00 





(B4) 


0 


v{x,y) y'^'^dy. 
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Here v{x,y) = C. 


a 


= 2\/r^ 


X = 


= \/(1^1 “ t)^ + s-nd Ks(a;) is the 


^^{x-R/2f+y^+r^ y/(x+R/2f+y^+r^ ^ 

MacDonald function. Unfortunately, function (IB4|) cannot be expressed in terms of elementary or special functions. 
However, it can be represented as a series 


Vs{x) 




dr = ^sign(a;) y^(— 1 )”+^ 

n—1 


(2n- 1)!! 
2"n! 


(B5) 


^s,n (- 2 ^) 


+ 00 

J e-"’’(r2 


1 


.^s-l/ 2 dr (2s- 1 )!!^ 3 o 

^ - 2«+i [T 


n+1/2 \ 

i, 0, (n - s) J 


(B6) 


where z = ax, a? = and Gf® is the Meijer G-function. When evaluating integral (IB6I) we used the formula 

2.1.4.20 from book and the relations of symmetry and shifting for the Meijer’s function. The condition that 
functions fk(x) and gk(x) are finite as x —>■ ±oo allows us to determine the spectrum. 
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